df <- read.csv("winequality-red.csv", header = TRUE, as.is = FALSE)
head(df)
## fixed.acidity volatile.acidity citric.acid residual.sugar chlorides
## 1 7.4 0.70 0.00 1.9 0.076
## 2 7.8 0.88 0.00 2.6 0.098
## 3 7.8 0.76 0.04 2.3 0.092
## 4 11.2 0.28 0.56 1.9 0.075
## 5 7.4 0.70 0.00 1.9 0.076
## 6 7.4 0.66 0.00 1.8 0.075
## free.sulfur.dioxide total.sulfur.dioxide density pH sulphates alcohol
## 1 11 34 0.9978 3.51 0.56 9.4
## 2 25 67 0.9968 3.20 0.68 9.8
## 3 15 54 0.9970 3.26 0.65 9.8
## 4 17 60 0.9980 3.16 0.58 9.8
## 5 11 34 0.9978 3.51 0.56 9.4
## 6 13 40 0.9978 3.51 0.56 9.4
## quality
## 1 5
## 2 5
## 3 5
## 4 6
## 5 5
## 6 5
Это данные о физико-химических свойствах вина, например, о кислотности, содержании диоксида серы, сульфатов и тд, а также качество вина (измеряемое баллами от 0 до 10).
1 - фиксированная кислотность
2 - летучая кислотность
3 - лимонная кислота
4 - остаточный сахар
5 - хлориды
6 - свободный диоксид серы
7 - общий диоксид серы
8 - плотность
9 - рН
10 - сульфаты
11 - алкоголь
summary(df)
## fixed.acidity volatile.acidity citric.acid residual.sugar
## Min. : 4.60 Min. :0.1200 Min. :0.000 Min. : 0.900
## 1st Qu.: 7.10 1st Qu.:0.3900 1st Qu.:0.090 1st Qu.: 1.900
## Median : 7.90 Median :0.5200 Median :0.260 Median : 2.200
## Mean : 8.32 Mean :0.5278 Mean :0.271 Mean : 2.539
## 3rd Qu.: 9.20 3rd Qu.:0.6400 3rd Qu.:0.420 3rd Qu.: 2.600
## Max. :15.90 Max. :1.5800 Max. :1.000 Max. :15.500
## chlorides free.sulfur.dioxide total.sulfur.dioxide density
## Min. :0.01200 Min. : 1.00 Min. : 6.00 Min. :0.9901
## 1st Qu.:0.07000 1st Qu.: 7.00 1st Qu.: 22.00 1st Qu.:0.9956
## Median :0.07900 Median :14.00 Median : 38.00 Median :0.9968
## Mean :0.08747 Mean :15.87 Mean : 46.47 Mean :0.9967
## 3rd Qu.:0.09000 3rd Qu.:21.00 3rd Qu.: 62.00 3rd Qu.:0.9978
## Max. :0.61100 Max. :72.00 Max. :289.00 Max. :1.0037
## pH sulphates alcohol quality
## Min. :2.740 Min. :0.3300 Min. : 8.40 Min. :3.000
## 1st Qu.:3.210 1st Qu.:0.5500 1st Qu.: 9.50 1st Qu.:5.000
## Median :3.310 Median :0.6200 Median :10.20 Median :6.000
## Mean :3.311 Mean :0.6581 Mean :10.42 Mean :5.636
## 3rd Qu.:3.400 3rd Qu.:0.7300 3rd Qu.:11.10 3rd Qu.:6.000
## Max. :4.010 Max. :2.0000 Max. :14.90 Max. :8.000
Удалим из рассмотрения free.sulfur.dioxide, так как есть total.sulfur.dioxide и удалим citric.acid, так как есть fixed.acidity.
df <- df[-3]
df <- df[-5]
Таблица внизу показывает, какую моду имеют наши количественные признаки. Отсюда можно сделать вывод о том, являются они непрерывными или нет. Построим matrix plot для того, чтобы увидеть особенности в наших данных.
library(lattice)
library(ggplot2)
library('GGally')
ggpairs(df, title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Прологарифмируем 1, 3, 4, 5, 8, 9.
dfl <- transform(df, fixed.acidity=log(fixed.acidity), residual.sugar=log(residual.sugar), chlorides=log(chlorides), total.sulfur.dioxide=log(total.sulfur.dioxide), sulphates=log(sulphates), alcohol=log(alcohol))
ggpairs(dfl, title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
Удалим единичные outliers.
dflo <- dfl
dflo[rownames(dflo)[dflo$volatile.acidity > 1.5 | dflo$chlorides < -4 | dflo$total.sulfur.dioxide > 5.1 | dflo$alcohol > 2.7],] <- NA
ggpairs(na.omit(dflo), title="correlogram", columns=c(1:9), upper = list(continuous = "points"), diag = list(continuous = "barDiag"))
names(dflo)[names(dflo) == 'fixed.acidity'] <- 'log_fixed.acidity'
names(dflo)[names(dflo) == 'residual.sugar'] <- 'log_residual.sugar'
names(dflo)[names(dflo) == 'chlorides'] <- 'log_chlorides'
names(dflo)[names(dflo) == 'total.sulfur.dioxide'] <- 'log_total.sulfur.dioxide'
names(dflo)[names(dflo) == 'sulphates'] <- 'log_sulphates'
names(dflo)[names(dflo) == 'alcohol'] <- 'log_alcohol'
Нормальность данных:
library(ggpubr)
ggqqplot(dflo$log_fixed.acidity, ylab = "log_fixed.acidity")
ggqqplot(dflo$volatile.acidity, ylab = "volatile.acidity")
ggqqplot(dflo$log_residual.sugar, ylab = "log_residual.sugar")
ggqqplot(dflo$log_chlorides, ylab = "log_chlorides")
ggqqplot(dflo$log_total.sulfur.dioxide, ylab = "log_total.sulfur.dioxide")
ggqqplot(dflo$density, ylab = "density")
ggqqplot(dflo$pH, ylab = "pH")
ggqqplot(dflo$log_sulphates, ylab = "log_sulphates")
ggqqplot(dflo$log_alcohol, ylab = "log_alcohol")
shapiro.test(dflo$log_fixed.acidity)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_fixed.acidity
## W = 0.98536, p-value = 1.234e-11
shapiro.test(dflo$volatile.acidity)
##
## Shapiro-Wilk normality test
##
## data: dflo$volatile.acidity
## W = 0.97928, p-value = 2.158e-14
shapiro.test(dflo$log_residual.sugar)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_residual.sugar
## W = 0.85681, p-value < 2.2e-16
shapiro.test(dflo$log_chlorides)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_chlorides
## W = 0.82868, p-value < 2.2e-16
shapiro.test(dflo$log_total.sulfur.dioxide)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_total.sulfur.dioxide
## W = 0.98794, p-value = 3.061e-10
shapiro.test(dflo$density)
##
## Shapiro-Wilk normality test
##
## data: dflo$density
## W = 0.99131, p-value = 4.189e-08
shapiro.test(dflo$pH)
##
## Shapiro-Wilk normality test
##
## data: dflo$pH
## W = 0.99327, p-value = 1.206e-06
shapiro.test(dflo$log_sulphates)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_sulphates
## W = 0.95749, p-value < 2.2e-16
shapiro.test(dflo$log_alcohol)
##
## Shapiro-Wilk normality test
##
## data: dflo$log_alcohol
## W = 0.946, p-value < 2.2e-16
Нормальность отсутсвует :( Но объем выборки позволяет говорить об асимптотической сходимости методов.
Посмотрим на корреляцию между признаками.
library(ggcorrplot)
library(ppcor)
dflo <- na.omit(dflo)
corr <- round(cor(dflo[c(1:9)]), 1)
ggcorrplot(corr, method = "circle", lab = TRUE)
ggcorrplot(pcor(dflo[c(1:9)])$estimate, method = "circle", lab = TRUE)
Сильная зависимость наблюдается между log_fixed.acidity и density, а также log_fixed.acidity и pH.
library(lm.beta)
lm.df.oecd <- lm.beta(lm(data = dflo, quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
summary(lm.df.oecd)
##
## Call:
## lm(formula = quality ~ log_fixed.acidity + volatile.acidity +
## log_residual.sugar + log_chlorides + log_total.sulfur.dioxide +
## density + pH + log_sulphates + log_alcohol, data = dflo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.67653 -0.35971 -0.05412 0.44655 1.91522
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 47.95847 0.00000 22.85206 2.099 0.036006 *
## log_fixed.acidity 0.49971 0.12368 0.21900 2.282 0.022635 *
## volatile.acidity -0.89892 -0.19787 0.10357 -8.680 < 2e-16 ***
## log_residual.sugar 0.10168 0.04479 0.06365 1.597 0.110370
## log_chlorides -0.22137 -0.08816 0.05872 -3.770 0.000169 ***
## log_total.sulfur.dioxide -0.07963 -0.06939 0.02490 -3.199 0.001408 **
## density -48.90313 -0.11373 23.17792 -2.110 0.035024 *
## pH -0.09064 -0.01738 0.19333 -0.469 0.639249
## log_sulphates 0.80965 0.22449 0.08466 9.563 < 2e-16 ***
## log_alcohol 2.63132 0.32207 0.29004 9.072 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6427 on 1582 degrees of freedom
## Multiple R-squared: 0.3639, Adjusted R-squared: 0.3602
## F-statistic: 100.5 on 9 and 1582 DF, p-value: < 2.2e-16
Значимыми являются с ***. согласно модели quality будет больше при меньшем значении 4 признаков (volatile.acidity, log_chlorides, log_total.sulfur.dioxide, density; максимальный по модулю– volatile.acidity) и большем значении 3 (log_sulphates, log_alcoho, log_fixed.acidityl; максимальный по модулю– log_alcohol). Остальные коэффицинты
Также выбор модели можно осуществлять на основании информационных критериев, которые базируются на функции максимального правдоподобия для оценок параметров и учитывают число параметров.
В случае функции stepAIC, AIC вычисляется по формуле \[\text{AIC} = 2k - 2 \log \mathcal L\], где \(k\) — число параметров модели, \(\mathcal L\) — функция максимального правдоподобия для модели.
“Backward” алгоритм stepAIC будет убирать переменные по одной, пока не будет достигнут минимум функции AIC.
lm.df.oecd.full <- lm.beta(lm(data = dflo, quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
stepAIC(lm.df.oecd.full, direction = "backward")
## Start: AIC=-1397.63
## quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar +
## log_chlorides + log_total.sulfur.dioxide + density + pH +
## log_sulphates + log_alcohol
##
## Df Sum of Sq RSS AIC
## - pH 1 0.091 653.55 -1399.4
## <none> 653.46 -1397.6
## - log_residual.sugar 1 1.054 654.51 -1397.1
## - density 1 1.839 655.30 -1395.2
## - log_fixed.acidity 1 2.151 655.61 -1394.4
## - log_total.sulfur.dioxide 1 4.226 657.68 -1389.4
## - log_chlorides 1 5.871 659.33 -1385.4
## - volatile.acidity 1 31.119 684.58 -1325.6
## - log_alcohol 1 33.998 687.45 -1318.9
## - log_sulphates 1 37.778 691.23 -1310.2
##
## Step: AIC=-1399.41
## quality ~ log_fixed.acidity + volatile.acidity + log_residual.sugar +
## log_chlorides + log_total.sulfur.dioxide + density + log_sulphates +
## log_alcohol
##
## Df Sum of Sq RSS AIC
## <none> 653.55 -1399.4
## - log_residual.sugar 1 1.449 655.00 -1397.9
## - density 1 3.576 657.12 -1392.7
## - log_total.sulfur.dioxide 1 4.144 657.69 -1391.3
## - log_chlorides 1 5.871 659.42 -1387.2
## - log_fixed.acidity 1 7.587 661.13 -1383.0
## - volatile.acidity 1 31.276 684.82 -1327.0
## - log_sulphates 1 38.914 692.46 -1309.3
## - log_alcohol 1 43.791 697.34 -1298.2
##
## Call:
## lm(formula = quality ~ log_fixed.acidity + volatile.acidity +
## log_residual.sugar + log_chlorides + log_total.sulfur.dioxide +
## density + log_sulphates + log_alcohol, data = dflo)
##
## Coefficients:
## (Intercept) log_fixed.acidity volatile.acidity
## 54.00449 0.58040 -0.90063
## log_residual.sugar log_chlorides log_total.sulfur.dioxide
## 0.11194 -0.21461 -0.07845
## density log_sulphates log_alcohol
## -55.26861 0.81480 2.56140
Информационный критерий AIC убирает всего 1 предикатор– pH.
Удалим признак log_fixed.acidity (большая корреляция со многими признаками) и построим модель линейной регрессии ещё раз.
corr <- round(cor(dflo[c(2:9)]), 1)
ggcorrplot(corr, method = "circle", lab = TRUE)
ggcorrplot(pcor(dflo[c(2:9)])$estimate, method = "circle", lab = TRUE)
Теперь не наблюдается сильных корреляций (как обычных, так и частных) между признаками.
lm.df.oecd.full <- lm.beta(lm(data = dflo[c(2:10)], quality ~ volatile.acidity + log_residual.sugar + log_chlorides + log_total.sulfur.dioxide + density + pH + log_sulphates + log_alcohol))
stepAIC(lm.df.oecd.full, direction = "backward")
## Start: AIC=-1394.4
## quality ~ volatile.acidity + log_residual.sugar + log_chlorides +
## log_total.sulfur.dioxide + density + pH + log_sulphates +
## log_alcohol
##
## Df Sum of Sq RSS AIC
## - density 1 0.057 655.66 -1396.3
## - log_residual.sugar 1 0.150 655.76 -1396.0
## <none> 655.61 -1394.4
## - log_total.sulfur.dioxide 1 5.390 661.00 -1383.4
## - pH 1 5.527 661.13 -1383.0
## - log_chlorides 1 7.075 662.68 -1379.3
## - volatile.acidity 1 34.832 690.44 -1314.0
## - log_sulphates 1 35.825 691.43 -1311.7
## - log_alcohol 1 72.703 728.31 -1229.0
##
## Step: AIC=-1396.26
## quality ~ volatile.acidity + log_residual.sugar + log_chlorides +
## log_total.sulfur.dioxide + pH + log_sulphates + log_alcohol
##
## Df Sum of Sq RSS AIC
## - log_residual.sugar 1 0.094 655.76 -1398.0
## <none> 655.66 -1396.3
## - log_total.sulfur.dioxide 1 5.335 661.00 -1385.4
## - pH 1 5.536 661.20 -1384.9
## - log_chlorides 1 7.147 662.81 -1381.0
## - volatile.acidity 1 34.812 690.48 -1315.9
## - log_sulphates 1 37.286 692.95 -1310.2
## - log_alcohol 1 108.388 764.05 -1154.7
##
## Step: AIC=-1398.03
## quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide +
## pH + log_sulphates + log_alcohol
##
## Df Sum of Sq RSS AIC
## <none> 655.76 -1398.0
## - log_total.sulfur.dioxide 1 5.250 661.01 -1387.3
## - pH 1 5.708 661.47 -1386.2
## - log_chlorides 1 7.056 662.81 -1383.0
## - volatile.acidity 1 34.720 690.48 -1317.9
## - log_sulphates 1 37.195 692.95 -1312.2
## - log_alcohol 1 112.828 768.59 -1147.3
##
## Call:
## lm(formula = quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide +
## pH + log_sulphates + log_alcohol, data = dflo[c(2:10)])
##
## Coefficients:
## (Intercept) volatile.acidity log_chlorides
## 0.36690 -0.93548 -0.23784
## log_total.sulfur.dioxide pH log_sulphates
## -0.08564 -0.43301 0.76832
## log_alcohol
## 3.10032
Информационный критерий AIC убирает 2 предикатора– density и log_residual.sugar.
Построим модель с использование минимального количества предикаторов.
lm.df.oecd <- lm.beta(lm(data = dflo[c(2:10)], quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide + pH + log_sulphates + log_alcohol))
summary(lm.df.oecd)
##
## Call:
## lm(formula = quality ~ volatile.acidity + log_chlorides + log_total.sulfur.dioxide +
## pH + log_sulphates + log_alcohol, data = dflo[c(2:10)])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.61118 -0.35630 -0.04643 0.45756 1.90911
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.36690 0.00000 0.50778 0.723 0.470065
## volatile.acidity -0.93548 -0.20592 0.10212 -9.161 < 2e-16 ***
## log_chlorides -0.23784 -0.09472 0.05759 -4.130 3.82e-05 ***
## log_total.sulfur.dioxide -0.08564 -0.07463 0.02404 -3.562 0.000379 ***
## pH -0.43301 -0.08301 0.11658 -3.714 0.000211 ***
## log_sulphates 0.76832 0.21303 0.08103 9.482 < 2e-16 ***
## log_alcohol 3.10032 0.37947 0.18774 16.514 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6432 on 1585 degrees of freedom
## Multiple R-squared: 0.3616, Adjusted R-squared: 0.3592
## F-statistic: 149.6 on 6 and 1585 DF, p-value: < 2.2e-16
Все коэффициенты регрессии значимы (***), согласно модели quality будет больше при меньшем значении первых 4 признаков (volatile.acidity, log_chlorides, log_total.sulfur.dioxide, pH; максимальный по модулю– volatile.acidity) и большем значении 2 оставшихся (log_sulphates, log_alcohol; максимальный по модулю– log_alcohol).
Анализ остатков:
library(MASS)
library(car)
par(mfrow = c(2, 2))
plot(lm.df.oecd, 2) #Normal Q-Q
plot(lm.df.oecd, 4) #Cook`s distance
r <- stdres(lm.df.oecd)
jackknife.res <- studres(lm.df.oecd) #Residuals vs Deleted
plot(jackknife.res ~ r, main = "Residuals vs Deleted") + abline(coef = c(0,1))
## integer(0)
Outliers не выявлено (значения Cook`s distance не превосходит 0.05, на графике “Residuals vs Deleted” все точки лежат на прямой y=x), остатки распределены нормально, значит построенная модель хорошо описывает зависимость.
Спрогнозируем значение quality для 2 индивидов с использованием построенной модели линейной регрессии. Первый индивид будет обладать значениями признаков близкими к средним, второй соответствует близким к “идеальным” значениям признаков (для данной модели).
new_wine <- data.frame(volatile.acidity = 0.52, log_chlorides = log(0.08), log_total.sulfur.dioxide = log(46), pH = 3.311, log_sulphates = log(0.65), log_alcohol = log(10.42))
new_wine2 <- data.frame(volatile.acidity = 0.015, log_chlorides = log(0.015), log_total.sulfur.dioxide = log(10), pH = 3, log_sulphates = log(2), log_alcohol = log(14))
predict.lm(lm.df.oecd, new_wine, interval = "confidence")
## fit lwr upr
## 1 5.654925 5.621168 5.688682
predict.lm(lm.df.oecd, new_wine, interval = "prediction")
## fit lwr upr
## 1 5.654925 4.392829 6.917021
predict.lm(lm.df.oecd, new_wine2, interval = "confidence")
## fit lwr upr
## 1 8.57 8.27884 8.861159
predict.lm(lm.df.oecd, new_wine2, interval = "prediction")
## fit lwr upr
## 1 8.57 7.275195 9.864805
Как и ожидалось первый индивид получил оценку quality близкую к средней по выборке. А второй получил оценку и вовсе больше максимального по выборке (уж сильно “идеальными” характеристиками он обладает).